NASA TN D-6066 


NASA TN D-6066 




r\ 




& 


COPY 






INVESTIGATION OF A DIGITAL AUTOMATIC 
AIRCRAFT LANDING SYSTEM IN XURBULENCE 


by Frank Neuman and John D. Foster 

Ames Research Center 
Moffett Field, Calif. 94035 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • OCTOBER 1970 





1 . Report No. 

NASA TN D-6066 


2. Government Accession No. 


3. Recipient's Catalog No. 


4. Title and Subtitle 

INVESTIGATION OF A DIGITAL AUTOMATIC AIRCRAFT LANDING 
SYSTEM IN TURBULENCE 


5. Report Date 

October 1970 


6, Performing Organization Code 


7. AuthorCs) 

Frank Neuman and John D. Foster 


8. Performing Organization Report No. 

A-3504 


9. Performing Organization Name and Address 

NASA Ames Research Center 
Moffett Field, Chlif. 94035 


10. Work Unit No. 

125-06-05-03-0021 


11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D. C. 20546 


13. Type of Report and Period Covered 
Technical Note 


14. Sponsoring Agency Code 


15. Supplementary Notes 


16. Abstract 


A digital system has been studied for automatically controlling the longitudinal motion of a large transport 
aircraft during the landing phase. The study was carried out by means of an all-digital simulation that was chiefly 
concerned with investigating the effects of gusts and wind shears on aircraft control near the ground. The 
performance of the autornatic control system operating in turbulence was determined by a Monte Carlo technique. 

With respect to the digital control system, it was found that (1) the basic analog flare mode could be modified 
to improve its performance under conditions of turbulence and wind shear; (2) for most control modes the 
computation rate requirement is surprisingly low, as indicated by the effects of computer repetition rate on the 
aircraft performance; (3) the performance degradation that results when a control computation cycle is occasionally 
skipped is relatively minor, a fact that is significant when the computer is shared with other systems for which 
emergency computations may have to be performed. 

Some of the causes of hard touchdowns showed the feasibility of developing a more efficient Monte Carlo 
technique that would aid future investigations. 


17. Key Words (Suggested by Author(s) ) 

Aircraft automatic landing 
Digital turbulence simulation 


18. Distribution Statement 

Unclassified —Unlimited 

19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

60 

$3.00 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 



















TABLE OF CONTENTS 


Page 

LIST OF SYMBOLS v 

SUMMARY 1 

INTRODUCTION 1 

SIMULATION OF THE AUTOMATIC LANDING SYSTEM 3 

Mode Controller 3 

Detailed Description of the Control Modes 5 

Altitude hold mode 5 

Capture mode 5 

Glide-slope tracking mode 6 

Flare mode 7 

REPRESENTATION OF ATMOSPHERIC TURBULENCE FOR 

AIRPLANE SIMULATION 10 

Background 10 

Simulation 13 

Digital simulation equations 14 

Properties of the turbulence model 16 

TECHNIQUES FOR ANALYZING THE SIMULATED FLIGHTS 17 

Statistical Analysis 17 

Paired Observations 19 

AUTOMATIC CONTROL SYSTEM PERFORMANCE 19 

Comparison of Analog and Digital Control System Performance 19 

Performance of the Noncritical Control Modes . . , 20 

Performance of the Critical Control Modes Under Steady Wind 

Conditions and Simple Shear 21 

Performance of the Critical Control Modes Under Turbulent Conditions 23 

Glide-slope tracking 23 

Rare mode 24 

Touchdown Performance of the Landing System . . . 29 

Performance of the System With Reduced Sampling Rates 30 

Performance of the System With Randomly Missing 
Computation Cycles 32 

CONCLUSIONS 33 

APPENDIX A - SIMULATION OF THE STABILIZED AIRPLANE 34 

AIRPLANE EQUATIONS 34 

Equations of Motion 34 

Engine and Engine Servos 39 

Elevator, Elevator Servo, and Raps 39 

iii 



TABLE OF CONTENTS 


Page 

AUTOMATIC FLIGHT-PATH STABILIZATION 40 

Pitch Control Loop 40 

Autothrottle 40 

APPENDIX B - SIMULATION OF LINEAR SYSTEMS WITH DIFFERENCE 
EQUATIONS 42 

APPENDIX C - FILTER GAIN CONSTANT ADJUSTMENT FOR THE DIGITAL 
TURBULENCE SIMULATION 48 


REFERENCES 


51 


LIST OF SYMBOLS 


drag coefficient, D/(l /2)pV^ S 
lift coefficient, L/(l/2)pV^S 

pitching-moment coefficient, pitching moment/(l/2)pV^Sc 

Cm at zero angle of attack 

wing mean aerodynamic chord, ft 

center of gravity 

aerodynamic drag, lb 

perpendicular distance of thrust line below center of gravity, ft 

X component of total aerodynamic force in the trajectory axis system, lb 
(appendix A) 

z component of total aerodynamic force in the trajectory axis system, lb 

filter transfer function 

acceleration due to gravity, 32.2 ft/sec^ 

reference sink rate for glide-slope extension, ft/sec 

altitude error from glide-slope beam center, ft 

initiation altitude of the feedback control law in the flare, ft 

desired vertical touchdown velocity, ft/sec 

sink rate at touchdown, ft/sec 

reference altitude for the flare law, ft 

altitude at which the predictive pitch ramp command begins, ft 
altitude at which the predictive pitch step command begins, ft 
instrument landing system, standard radio guidance installed at major airports 
aircraft moment of inertia around y-axis, slug-ft^ 



ij 

k 

ke 

L 

Lu,Lw 

m 

n 

"h 


nx 


thrust angle of incidence with body respect to x-axis (up, positive), radians 

filter gain constants 

glide-slope tracking gain constant 

aerodynamic lift force, lb 

horizontal and vertical turbulence scale lengths, respectively, ft 
airplane mass, slugs 

number of observations in a statistical sample 

percentage of li feedback as compared to h feedback in the flare law 
x-component of gust in earth axis 


nz 


q 

Ru(r)jRw(r) 


S 


s 

s 


2 


T 

t 

V 

Vc 


Vx 

Vx 

Vz„ 


g 


g 


Vz„ 


z-component of gust in earth axis 
pitch angular velocity, radians/sec 

spatial autocorrelation function of the longitudinal and vertical components of 
turbulence, respectively 

wing reference area, ft^ 

Laplace transform operator 

sample variance 

total thrust, lb 

sampling time interval, sec 

airspeed, ft/sec unless otherwise indicated 

velocity command, ft/sec 

velocity of the aircraft with respect to the inertial frame of reference, ft/sec 
x-component of total wind in earth axis, ft/sec 
x-component of mean wind in the earth axis system, ft/sec 
z-component of total wind in earth axis, ft/sec 
z-component of mean wind in the earth axis system, ft/sec 


VI 



X 


o 


z 

Of 

7 

^ILS 

5e 

5f 

e 

V 

e 

M 

p 

o 

n 

CO 


horizontal distance of aircraft at touchdown from the ground intercept of the ILS 
glide-slope beam (touchdown distance) 

z transform operator 

ground effect 

angle of attack, radians 

flight-path angle, radians 

ILS glide path angle, radians 

elevator deflection, radians 

flap deflection, radians 

glide-slope error angle, radians 

capture initiation angle, radians 

angle between body axis and trajectory axis of the aircraft 
pitch angle of airplane body axis relative to horizon, radians 
sample mean 
air density, slugs/ft^ 

standard deviation of the quantity denoted by the subscript 
spatial and temporal power spectral density of the turbulence 
spatial frequency of the gusts, radians/ft 
angular frequency, radians/sec 


Subscripts 

c command signal 

d desired quantity 

i input, or initialization quantity 


vii 


o 


output, or quantity at touchdown 


component tangential to the flight path 
longitudinal and vertical components of the quantity 
average 

first and second runs of a paired observation 


INVESTIGATION OF A DIGITAL AUTOMATIC AIRCRAFT LANDING SYSTEM 


IN TURBULENCE 
Frank Neuman and John D. Foster 
Ames Research Center 


SUMMARY 


A digital system has been studied for automatically controlling the longitudinal motion of a 
large transport aircraft during the landing phase. The study was carried out by means of an 
all-digital simulation that was chiefly concerned with investigating the effects of gusts and wind 
shears on aircraft control near the ground. The performance of the automatic control system 
operating in turbulence was determined by a Monte Carlo technique. 

With respect to the digital control system, it was found that (l)the basic analog flare mode 
could be modified to improve its performance under conditions of turbulence and wind shear; 
(2) for most control modes the computation rate requirement is surprisingly low, as indicated by 
the effects of computer repetition rate on the aircraft performance; (3) the performance 
degradation that results when a control computation cycle is occasionally skipped is relatively 
minor, a fact that is significant when the computer is shared with other systems for which 
emergency computations may have to be performed. 

Some of the causes of hard touchdowns showed the feasibility of developing a more efficient 
Monte Carlo technique that would aid future investigations. 


INTRODUCTION 


Recent history of aviation has seen a continuing trend toward automation. Most current 
transport aircraft are equipped with analog automatic control systems that are used routinely during 
the cruise phase and can be coupled to the ILS beam for portions of the approach. Systems have 
now been developed that will permit automatic control of aircraft to touchdown. Although research 
aircraft use these systems for landing under zero visibility conditions, for commercial operation 
there is still concern for the accuracy of the guidance information, the reliability of the automatic 
systems, and the effectiveness of situation information displays necessary for recovery from 
malfunctions. 

This investigation was undertaken because it was believed that digital control techniques could 
contribute to the objective of achieving acceptable automatic landings under all operational 
conditions. Four specific contributions were identified: (1) Combining many functions in a 
computer relieves the pilot from secondary manual tasks now associated with automatic landing 
operations. The ultimate goal would be for the pilot to become a decision-making systems’ manager 
rather than a primary control element. (2) Data stored in the computer are available for selective or 



automatic call and observation by the pilot. This flexibility is also important in a research 
environment where one attempts to define the type and form of the data needed by the pilot in his 
new role as systems’ manager. (3) Digital computations are repeatable precisely. Since redundant 
systems are required for safety, several computers making calculations based on identical sensor 
data must produce identical answers. Failure to do so is an instant warning of malfunction. (4) More 
sophisticated control laws can be implemented that will improve performance at a small cost of 
computer storage and computation time. 

While digital control has been used extensively in military aircraft, the economic factor has so 
far prevented its commercial application. However, since transport aircraft projected for the future 
will have a substantial on-board digital computer capability for other purposes for which analog 
computations are unsatisfactory (such as navigation computations, fuel management, and engine 
performance monitoring), economics will no longer be a deterrent. 

Flight control requirements for commercial jet transports differ significantly from those of 
military aircraft. It is therefore appropriate to examine digital control problems specifically for 
commercial jet transport aircraft. A primary requirement is extreme safety. One aspect of flight 
safety is system reliability. Another aspect is the effect of the environment such as noisy guidance 
information and wind turbulence. System reliability has been studied extensively and various 
techniques, such as component redundancy and self-checking, have been investigated in depth. It 
has been shown in reference 1 that turbulence causes larger deviations from the desired flight path 
than the errors in ILS guidance. This study therefore concentrated on the effect of turbulence on 
safe automatic landings. 

The present study had three specific objectives. The first was to choose an efficient method of 
translating analog autopilot technology directly into digital control laws to take advantage of 
previous experience in analog autopilot design. This was a desirable step before considering 
improved designs which take advantage of the superior logical capability of the digital machine. The 
second objective was to study some improvements of an analog flare law. This goal was easily 
accomplished in the digital system but would have been difficult to implement in an analog system. 
The third, and main objective of the study, was to determine the influence of wind turbulence on 
the performance of the automatic control system during landing. The very nonlinear problem of 
flare and touchdown was attacked by Monte Carlo methods. If primary attention had been focused 
on the approach phase, which can be effectively expressed by a linear model, the more efficient 
power spectral method with gaussian noise, described in reference 1, could have been used. 

In line with the above objectives, a digital method was developed to generate the turbulent 
wind components. Since the method has wide application, it is described in detail in the body of the 
report. 

For this investigation, the scope of the automatic landing problem was restricted in two ways. 
First, the jet transport simulation equations reported in reference 2 were reduced to three degrees 
of freedom by considering the longitudinal axis only. This restriction is reasonable in light of the 
accident statistics compiled in reference 1 which concludes that accidents due to longitudinal errors 
are fatal much more often than accidents due to lateral errors. Second, the system guidance 
information was assumed to come from an error free ILS beam and altimeter. 
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SIMULATION OF THE AUTOMATIC LANDING SYSTEM 


In this section the longitudinal automatic landing system will be described, and some of the 
design considerations will be given. The system is illustrated by combinations of flow charts and 
block diagrams that indicate its digital nature and its analog heritage. Reference 3 is an excellent 
discussion of linear analog aircraft control laws. Additional information of modern autopilot design 
was obtained by private communication from the author of reference 4. While linear filters will be 
described by transfer functions, it should be understood that the filtering function is actually 
performed in the digital computer by solving difference equations. This process is described in 
appendix B. 


Mode Controller 

The flight control laws are segmented into control modes for different portions of the 
approach and landing. A mode controller automatically selects the proper control mode in sequence 
(i.e., the altitude hold, glide-slope capture, glide-slope tracking, and flare mode) according to 
predetermined criteria. In accordance with conventional practice (ref. 4), the control modes operate 
on the velocity and pitch-stabilized aircraft (see fig. 1 ), and therefore operate with only three 
command variables, velocity command V^, pitch angle command Oc, and flap command 5f^. Data 
on the basic airplane and design of the inner stabilization loop are contained in appendix A. 


I 1 



L 


I 

Digital airplane simulation 
(pitch stabilized} 


Figure 1 .— Overall block diagram of the simulation. 
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The operation of the mode controller is shown in the form of a Fortran flow chart in figure 2. 
While the controller selects from only four modes, it passes through six different conditions. 
Conditions 1 and 4 are preparatory to actual mode switching. 


Go to {1 , 2, 3, 4, 5, 6) Cond 
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Figure 2 — Mode control logic. 


Altitude Glide slope Glide slope 



Figure Automatic landing geometry using ILS and 
radar altimeter. 


The mode controller is best described by 
considering a landing approach (see fig. 3). 
The aircraft approaches the ILS glide slope at 
a constant altitude of 1500 feet (condition 1, 
fig. 3). It penetrates the beam until the 
glide-slope error angle detected by the ILS 
receiver is equal to or less than 0.5^. At this 
point the capture initiation angle ec is 
calculated from the aircraft altitude h and 
the glide-slope angle Tjls- 

"c = ^ILs[»-«''0’*2500 

The mode controller then goes to condition 2 
until ec is reached. This permits the 
initiation of capture at the same distance 
from the beam center (2500 ft) independent 
of aircraft altitude. At this point the flaps are 
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commanded to extend to 50° and the controller switches to the capture mode. The capture mode, 
condition 3, is timed and switched to glide-slope tracking, condition 4, after 5 seconds. Glide-slope 
tracking proceeds to a preselected altitude, at which point the sink rate and velocity of the aircraft 
are used to calculate the flare initiation altitude h^j^p and other initiation parameters described in 
the flare mode section. The mode controller is then switched to condition 5 and glide-slope tracking 
is continued until h^j^gp is reached, at which point the flare mode is selected. 

Notice that no automatic go-around mode is provided. The simulated aircraft is forced to land 
so that the conditions can be found that result in unsatisfactory landings. 


Detailed Description of the Control Modes 

Altitude hold mode— A simple 
altitude hold mode incorporated in the 
system provides the initial conditions for 
capture. The digital control was modeled 
after the representative analog system 
shown in figure 4. The system consists of a 
differencing circuit for calculating the 
altitude error h, followed in series by a low 
pass filter, a gain, and a low gain integrator. 
For comparison, the digital equivalent 
equations are as follows (see eq. (8) 
table 8): 


0c = + C20^_^ + C3Ah + C4Ah_^ 

These equations are solved once each computation cycle. 



Typical value: k - 0.00128 ft’’ 8 ^ 

Figure 4,- Altitude hold mode. 


Ah = h - h ^ 
ref 



Typical value; k = -0.0005 sec/ft 8 

^Cmin = -11.4 ft/sec 


Figure 5.- Capture mode. 


Capture mode— The capture mode 
(fig. 5) provides for a smooth rotation from 
level flight to the glide-slope angle. This 
mode, specifically designed for this digital 
simulation, differs from capture modes 
used in analog systems in that a small step 
pitch angle command A0p is applied at 
capture initiation to rotate the airplane. 
The magnitude of the step is based on the 
glide-slope angle of the beam to be 
captured. In addition, an inertial vertical 
velocity error signal is generated to increase 
the sink rate linearly from capture 
initiation to the proper sink rate for the 
given glide-slope angle and velocity 
command. The pitch error signal 0^ is 
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then integrated and filtered to produce the pitch angle command. The integrator provides an error 
signal proportional to altitude error. Since the sink rate reference is a ramp function whose final 
value is the proper sink rate of the aircraft on the glide slope, the resulting altitude reference is a 
parabolic curve that smoothly intersects the glide slope. As shown in the mode control logic at the 
beginning of capture the flaps are deployed from 35° to 50° causing a smooth reduction of 

airspeed. 

Glide-slope tracking mode- After glide-slope capture the aircraft remains in the glide-slope 
tracking mode until flare. The glide-slope error angle e is filtered as shown in figure 6. The tracking 
loop operates on both the error and the integral of the error. Its gain is adjusted proportional to the 
altitude from capture to 200 feet. Reducing the gain with decreasing altitude will correct signal 
output approximately in proportion to displacement of aircraft from the beam center instead of in 
proportion to angular displacement measured by e. Because of various signal reflections in the 
region close to the ground, the path becomes so irregular that it is difficult to follow. Therefore, the 
above gain reduction is insufficient and kg is smoothly reduced to zero between 200 and 1 00 feet 
altitude. In air turbulence the aircraft tends to wander around the glide-slope center, occasionally 
experiencing large sink rates. These are a problem only at altitudes below 200 feet where it is 
important to keep the sink rate close to its proper value so that the flare mode can safely land ^e 
aircraft. Therefore, a glide-slope extension signal added to the tracking signal damps the tracking 
response below 200 feet and eventually controls the aircraft to a constant sink rate regardless of its 
glide-slope position. This extension is a pitch attitude command Oc, proportional to H - hj,, 
where h|j is the proper sink rate for the given glide slope. 


Glide slope tracking 



fib = -Vc y 


ILS 


kv = -0.004 


Glide slope extension 


Figure 6.- Glide-slope tracking mode. 
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Flare mode- The flare mode controls the traditional exponential flare (ref. 4). The flare has 
three boundary conditions: 


hj- initiation altitude of the feedback control law 
fif initial sink rate 

^fo desired vertical touchdown velocity (taken as 1.5 ft/sec for this simulation) 


A flare law that satisfies these boundary conditions is: 


-t/a2 

hj.(t) = (hf - a2h^^)e + a2h^^ 


and a 2 is calculated as 


^2 = -hf/(h^ - 

The reference sink rate is the derivative of equation (1) 




CD 

(la) 


-t/a. 


( 2 ) 


hr = -(l/a2)(h^ - a2h^^)e 

The predictive portion of the flare law (fig. 7) has two sections, a step command in pitch, A0p, 
which causes the aircraft to begin to rotate, and a ramp pitch command A0j^, which begins 
somewhat later. With no other disturbance, the predictive flare commands will generate an 
approximately exponential flare. Feedback is used to overcome disturbances. Equation (1) is the 
solution of the following differential equation: 

hr + ^2Chr + hfo) (3) 


Corrective 

feedback 

signal 

generation 



Typical values; ki = -0.0012 
ai = 0.333 
32 = 5.4 

Figure 7.— Flare mode block diagram. 
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with the boundary conditions hr = hf at t = 0 and hr - hfo at h - 0. The feedback portion of the 
flare law generates a corrective signal when equation (3) is not fulfilled by the actual altitude h and 
sink rate h in place of hr and fir- The corrective signal is (see fig. 7). 

6^ = k£[l + (ai/s)][h + a 2 (h - h£)] (4) 

which is added to the predictive pitch command. Hence, no correction signal is applied when the 
reference path is followed. 

There is, however, a difference between the flare laws implemented on modern analog 
automatic landing systems and the flare law implemented here. In many modern analog systems the 
predictive pitch commands as well as the feedback gain constants and the flare initiation 
altitude hr are preselected for the nominal glide-slope angle. In our digital system these values are 
computed from the flight-path angle 7 just before flare initiation. This method stabilizes the 
landing performance when the flight path is disturbed by windshear or turbulence. 

The flare law used in this investigation is shown in detail in the form of a Fortran flow chart in 
figure 8. When the flare subroutine is entered for the first time, the sink rate is used to calculate 
decision altitudes for the predictive flare law commands. The altitude at which these initial 
calculations are made is somewhat above the highest at which the flare may be started. As figure 8 



Figure 8.- Fortran flow chart of flare computations. 
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shows, the step command altitude is proportional to the flight-path angle y ^ -li/V. The ramp 
begins at a proportionally lower altitude. Since the aircraft docs not begin to deviate from a 
straight-line glide path instantaneously upon receiving the pitch step command, the altitude for the 
corrective feedback to begin is also selected proportionally lower than the step command altitude. 
After these calculations are completed, the flare computer transfers the authority back to the 
glide-slope tracking control. 

When the step command altitude, h^tep. is reached, the flare control mode takes over 
completely. First a thrust reduction ATc is calculated and from the latest measured sink rate, h, 
the step command A0p is calculated and executed. Then the ramp increment and the 

feedback gain constant ^2 calculated. At this point the flare computer is switched to its final 
submode. 

In the final submode the predictive ramp pitch command is added to the corrective feedback 
flare command'. The summed signals are transmitted as the pitch change command 0c to the pitch 
control loop of the aircraft. In addition, a timed thrust reduction is programmed to reduce the 
speed at touchdown. 

Some features of the flare mode of figure 6 were incorporated because of the simulation and 
flight-test results of a flare system reported in reference 5. In that system (1 ) command rate limiting 
was used, (2) the corrective feedback integrator gain was positive for pitch-up commands only and 
zero for pitch-down commands, and (3) the total pitch command could never be smaller than the 
initial pitch command. The first two features added to the present flare law improved the 
performance by reducing touchdown velocities and thus were incorporated in the final simulation. 
The best rate limiter was found to be 0.0926 radian per second. The best integrator gains were 
achieved when the pitch-down gain was half the pitch-up gain. The third feature did not improve 
performance so was not used. 

Under disturbances, the feedback term in the flare law (eq. (4)) does not attempt to guide 
along a path fixed in space, or even hold h(t) and h(t) at given values. As long as the feedback 
signal of equation (4) is zero no correction is made. Disturbances, therefore, tend to cause 
translations of the touchdown point rather than large maneuvers to meet a given touchdown point 
which would often cause hard landings. 


With the above remarks in mind, an improvement in the predictive flare command was made, 
which is not shown in figures 6 and 7. Instead of a ti med ramp signal, the total increase in pitch 
command due to the ramp was applied between the ramp initiation altitude (hi-^j^ip) and ground as 
a function of change in altitude between computation cycles, hn - hn-i- Thus the change in pitch 
command for each cycle is computed as: 


A0 - 


This change replaces the timed ramp increments. Thus, when 
altitude, the predictive pitch command is applied at 
change remains constant. 


wind shear causes a fast drop in 
a faster rate, while the total 
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Another measurable variable h often available from aircraft sensors might improve the flare 
performance of the airplane. A differential equation that contains h and that also has equation (1 ) 
as a solution is 

h + n^a2h - (1 - n^)a2^h - a2h£ = 0 (5) 

where nj^ is a factor between 0 and 1, and nj^ = 1 results in the previously described feedback 
signal of equation (4). If vertical acceleration is measured, equation (5) with properly 
chosen nj^ may provide the corrective signal that will improve performance under turbulent 
conditions. 


REPRESENTATION OF ATMOSPHERIC TURBULENCE FOR AIRPLANE SIMULATION 


Background 

Before presenting the analysis and results of this study, it is important that the atmospheric 
turbulence model be understood, because it has a major effect on the experimental results. The 
following section therefore presents the mathematical development of the digital turbulence model. 

Turbulence has three velocity components. For the pitch axis control problem only the 
longitudinal and vertical gust components ug and wg are of concern. Both are functions of both 
time and position. Figure 9 shows that the measured turbulence records obtained from an aircraft 
will be different for different flight speeds over the same flight path. Therefore, the spatial spectra 
calculated from the flight records will be different. Also, because of the limited sample time, the 
measured spectra differ from the actual spectra. For reasonable flight speeds, the changes in ug and 
Wg are smaller with respect to time than with respect to position. Therefore, for purposes of 
calculation, the gust velocity vectors at each point will be assumed not to vary with time. This is the 
concept of the turbulence field frozen with respect to time. Then, as figure 10 shows, the time 
records for two flights at different but constant speeds will be similar to each other and to the 





Figure 9 Time and position history of vertical turbulence. 
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Fast 

flight 



Figure 10.— Frozen gust pattern of cog. 

spatial distribution. If the flight speed is known at all times, even if it is varying, the time record can 
be converted to a spatial record of turbulent velocity components. To obtain the parameters of the 
mathematical model of turbulence, the spatial records are then fitted to the theoretical power 
spectra. 


A particular set of power spectra, the Dryden spectra (ref. 6), is easy to simulate and shows 
reasonable agreement with measured spectra and was therefore selected for this simulation. For the 
Dryden spectra one chooses the longitudinal spatial autocorrelation function of gust velocity to be a 
simple exponential (see ref. 5) 




-r/Lu 


( 6 ) 


From the equation of continuity, Batchelor (ref. 7) has shown that the longitudinal and vertical 
autocorrelation functions are related by the equation 






r 9Ru(r) 
2 3r 


(7) 


Therefore, since for homogeneous turbulence Oy - Lu - L^, 



( 8 ) 


Tlie one-dimensional spatial power density spectra are obtained by substitution 
(8) into the well-known relation 




2 

TT 



R(r)cos (nr) dr 


of equations (6) and 


C9) 


1 1 



resulting in 


^ 


1 + 


= O 


2 ^ 1 + 3Lw^!^^ 
^ ^ (l + 


( 10 ] 


In these and the following equations 
L scale length, ft 

i2 spatial frequency, rad/ft 

CO angular frequency, rad/sec 

V aircraft ground speed, ft/sec 

a root-mean-square value of turbulence velocity, ft/sec, standard deviation 

r correlation interval, ft 

It is also assumed that the horizontal and vertical power spectra are uncorrelatcd. 

The actual turbulence depends on the two parameters Land a. Insufficient data exist to 
define L and o for all conditions. In the interim the following relationships will be used. 


Ou = 0.2 Ug 

= 0.2 Ug(0.5 + 0.00098 h) 0 < h < 500 ft ► 
= 0.2 Ug h > 500 ft , 


( 11 ) 


where Ug is the mean longitudinal wind velocity and h is the altitude above ground. Also the 
horizontal gust scale length is much less affected by altitude than is the vertical scale length. 


Lu = 100 Vh h > 230 ft ' 
Lu =600 h < 230 ft 

Lw = h ' 


(11a) 


In general, turbulence increases with mean wind. Because of the effect of the ground there is 
also an increase of vertical gust component with altitude (ref. 8, fig. 4) while the mean wind and 


r 
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horizontal gust component are independent of altitude. These relationships apply for approach 
altitudes up to 2000 feet, which are considered in this report. Thus, on the ground, wing altitude is 
10 feet, Lu is 600 feet, and Lw is 10 feet, but they are equal at 1000 feet. 

It should be noted at this point that the relationships for Lu, Lw, uu. and are not 
consistent with the assumption of homogeneous turbulence. However, it is assumed throughout the 
literature that using homogeneous equations with different values for the horizontal and vertical 
parameters provides a more realistic description of turbulence. 


Simulation 


For simulation purposes the spatial spectra are converted to temporal frequency spectra, 
which are a function of V, the aircraft velocity (ref. 5), 

O(u)) = i (12) 

Using equation (1 2), and noting that L2 = co/V, equations (10) become 


5’uCw) 



2^^ 1 

1 + (Lu/V)^w^ 


( 13 ) 


<I>w (w) 



1 + 3(Lw/V)^g)^ 

"" [l + (Lw/V)2o)2]2 


( 14 ) 


The above equations are defined so that 

o 

for equations (13) and (14). Note that the frequency spectra throughout this section are defined for 


For simulation it is necessary to generate turbulence with the 
required spectrum 4>(co) of standard deviation o and scale 
length L for a given flight velocity and altitude. A wide band 
noise thus applied to a linear filter of the proper 

frequency response to obtain 4>y(co) (see fig. 11). In order to carry 
out this operation, first find the filter transfer function F(s). The 
relation of input to output power spectrum is 

^y(u)) = |F(s)|^^.^$x(<^) (15) 

If the input power spectrum is white noise <I>x(co) = 1 , then 

^y(o)) = |F(s)|J_ (16) 

' I I S-J 0) 


positive frequencies only. 



Figure 1 1 .— Power spectrum 
shaping. 
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Now define 


V 

- Lu ’ 


/3 Lw 


ku = 


kw “ 




3awOw‘^ 


In conjunction with equation (16), equations (13) and (14) become 

k 


|Fu(jw) 


2 _ 




+ au^ 


iFwCjt^) 


which results in 


Fu(jw) = 


ju + au 


FwCjw) = 


2 _ kw(u)^ + b^) 
(u)^ + aw^)^ 


Cj“ + b) 


(jo) + a„)' 


(17) 


Finally for 



(18a) 

F„(s) - ^ -4^ 

(18b) 

(s + aw)^ 



Digital simulation equations^ Since the noise is to be generated digitally, one must supply a 
digital equivalent of random gaussian noise and either write the differential equations from the 
trSisfer functions and solve them by numerical integration, or develop difference equations 
approximate the filters. The filters of equations (18a) and (18b) were simulated by means o 
dSerence equations - since this is an efficient method of computation. The principles used to 
tnerate the difference equations are discussed in appendix B. From table 8 and equation (18b), the 
difference equation for the vertical spectrum is 

Tn = CiXn-i - C27n-2 + DiXn-i + D2Xn-2 U9) 

where the constants of equation (19) are 


Cl = 2e 




Co — - e 


-2awT 


( 20 ) 

( 21 ) 



( 22 ) 


( 23 ) 
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( 24 ) 


The difference equation for the horizontal spectrum described by equation (1 8a) is 

Xn = Cyn-i + Dxn_i 
where the constants of equation (24) are 


C = 


e 


-au 


T 




(25) 

(26) 



rs 


Figure 12.— Digital simulation of turbulent gust. 


The filter input is a gaussian random 
sequence (see fig. 1 2). A good approximation 
of the wide band gaussian noise input 
sequence is provided by summing five pseudo 
random numbers (rj) with a flat distribution 
in the range -1 to +1 (see, e.g., ref. 9) 


zi 




£ 

3 


and adjusting the variance of the sequence to generate the proper output variance by selecting the 
multiplication factor k (xj = kzj, For this purpose the ratios of actual output to 

input variance ^^yg/^^xg the filters in figure 1 2 must be found. A good approximation for 
small T has been calculated for both filters in appendix C: 


a _ 


T 2 

IT y 


(27) 


where Oy^ is the desired and Oy^ the actual output variance of equation (13) or (14). Hence, to 
generate the desired output variance, Oy^ = 0 ^ 3 ^, the variance of the input sequence, ax^ must 
be = tt/T. This is achieved by premultiplying each Zj by the factor k to give Xj, 


k 



(28) 


For efficient calculations the constant k of equation (28) is combined with the input 
coefficients Dj of the difference equations (19) and (24). 
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Figure 13.- Reduction of output variance of the vertical 
power spectrum with increase in sample time 
interval (V = 220 ft/sec). 



Figure 14.- Autocorrelation of the vertical gust 
spectra. 


Properties of the turbulence 
morfd- Figure 13 shows the reduction in 
output variance due to increasing T as 
calculated from the exact equation in 
appendix C. The reduction is greatest for 
short scale lengths (low altitudes) and long 
sampling intervals because the 
high-frequency components are not 
presented. However, since the rigid aircraft 
does not respond to high frequencies, the 
absence of them does not affect the results 
of the simulation. 

To verify that the digital turbulence 
simulation behaved as predicted from the 
Dryden spectra model, 4000 noise samples 
for different scale lengths (L= 200, 1000, 
5000) for t = 0. 1 sec were generated and 
the autocorrelations calculated. The 
autocorrelation functions of the limited 
samples agreed well with the theoretical 
prediction (see fig. 14). The adequacy of 
generating noise samples at 0. 1 sec intervals 
was verified by statistical runs made with 
the simulated aircraft at t = 0.025 sec and 
at t = 0. 1 sec, the value normally used in 
this simulation, and no significant 
difference in systems performance was 
found. 

An important property of the 
turbulence is the change in turbulent 
velocity, Aug, over a given 
distance r because such changes affect the 
flight path of the airplane. This gradient 
property for the longitudinal Dryden 
spectrum has been derived in reference 10 
for a random initial value. Since Ug is 
gaussian, so is Aug by symmetry with a 
mean of zero. The variance of Aug is 


2 



and from equation (6) 


[Ug(x - r) - Ug(r)] = 2[ou^ - Ru(r)] 


(29) 


2 

a, 

AUg 



e 


-r/Lu 


) 


(30) 
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Equations (29) and (8) can be used to derive the gradient properties of the vertical gust component, 

V, " (31) 

o 

The gradient properties of the simulated horizontal turbulence were checked for a 2,000,000-foot 
sample with airspeed and scale length typical of landing and a variance of 16. The cumulative 
distributions Aug at various separations r for Lu = 100 to Ly = 1000 feet were plotted on normal 
probability graph paper (not shown) and the measured gradients agreed well with the theoretically 
predicted gradients. 

In summary, for the digital simulation, equations (1 1) and (1 la) were applied to calculate the 
turbulence parameters as functions of altitude and average wind velocity. These parameters were 
then taken to calculate the difference equation coefficients (eqs. (20)-(23), (25), (26)) whenever the 
altitude of the aircraft changed. Finally, a gaussian random sequence of Xg’s (fig. 1 2) of the proper 
variance (eq. (27)) was generated and inserted into the difference equations (19) and (24) to 
calculate the vertical and longitudinal gust components, y^. 


TECHNIQUES FOR ANALYZING THE SIMULATED FLIGHTS 


Statistical Analysis 

Under turbulent conditions, system performance has to be measured statistically. While it is 
possible to measure errors for each control mode, performance indexes were established only for 
the glide-slope tracking mode and the flare mode, because they were the major influence on the 
automatic landing performance. Information given in the section on results indicates that the 
altitude hold and capture modes performed satisfactorily. 

The glide-slope tracking mode provides the proper initial conditions for the flare mode. To 
track correctly, the airplane must remain well within the beam, except during glide-slope extension. 
The ability to remain within the beam was measured by two criteria. The first measured the mean 
and standard deviation of the altitude error from the center of the ILS beam, hg, These values 
were calculated for each flight from the end of capture to the beginning of the flare. The altitude 
error instead of the glide-slope angle error was chosen as a guidance measure even though the 
guidance signal is proportional to the glide-slope angle because the other choice would have 
weighted path errors close to flare very heavily. The second criterion, which also measured the 
ability to remain within the ILS glide-slope beam, was taken from an FAA advisory circular 
(ref. 1 1 ), which deals with automatic-pilot-coupler systems used as part of a category II installation. 
It requires system performance criteria to be met under the following wind conditions: a 
downwind component of 10 knots and, commencing at an altitude of 500 feet, a wind shear of 
4 knots per 100 feet altitude. To provide satisfactory glide-slope performance the criterion requires 
that the aircraft be stabilized on the glide slope before reaching an altitude of 700 feet above the 
field level. From 700 feet altitude to the decision altitude the autopilot-coupler should cause the 
airplane to track the center of the indicated glide slope to within ±35 microamperes or ±12 feet, 
whichever is the larger, without sustained oscillations. Notice that the criterion compensates for the 
increased position error sensitivity at low altitudes by changing from an angular error 
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(1 50 /lA = 0.7°) to an altitude error. While the criterion is apparently to be applied only for the 
windshear condition stated, it was applied to all simulated flights in this investigation, and the 
percentage outside the above limits was calculated. 

For a successful flare, reasonable initial conditions must be established. As a measure of these 
conditions, which are the terminal conditions of the glide-slope tracking, the airplane flight variables 
were recorded at 75 feet altitude, and for each series of flights the mean and variance of the 
parameters were calculated. Only the significant parameters for a successful flare will be discussed. 

For the flare performance only the sink rate at touchdown, ho, and touchdown distance from 
the glide slope to ground intersection, Xq, are calculated. Since we are not actually controlling for 
minimum deviation from a desired flight path, mean squared flight-path error was not calculated for 
the flare. 


Landing point dispersions caused by the flare control law can be separated from those caused 
by horizontal position errors at flare initiation. Since the flare law does not operate on horizontal 
position errors, the variance of the landing distance from glide slope ground intercept is the sum of 
the flare law induced errors and the flare window position errors. We define the flare window 
position error as the aircraft position error from the glide-slope position at 75 feet altitude. 
Computer time is saved because it is possible to initialize the aircraft on the glide slope at an 
altitude of 200 feet when one intends to study effects of altering the flare law. Starting the system 
from 45 000 feet out results in the flare window position error, which is characteristic of the overall 
landing system. Starting at an altitude of 200 feet on the glide slope (3800ft from glide slope 
intercept) results in a smaller flare window position error. If the smaller flare window position error 
is subtracted from the touchdown position error and the characteristic flare window position error 
for the overall system is added, the proper overall touchdown dispersion is 


|xi=45,000 


2 

= 


Xi=3 , 800 


h=75 

Xi=45,000 


2 

^X 


h=7 5 
Xi=3,800i 


The effects of altitude errors were insignificant, since starting the aircraft simulation from 200 ft 
altitude introduced sufficient attitude dispersion at the 75-foot altitude. An example is given using 
data for the case of 20 feet head wind. 

= (540)2 + [(120)2 _ (31^2] = (546)2 
As one can sec, the effect of neglecting the second term would be rather small. 

For flights without turbulence only one flight was necessary for each condition since the 
digital computer repeated identical flights exactly. In turbulence at least 100 flights were simulated 
for each of seven different values of gust variance. For each series of flights the mean and standard 
deviations of all the performance measures were calculated (e.g., he, ohg). 
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Paired Observations 


When a comparison is made between two slightly different control laws in an actual noisy 
environment, the variation of the noise may mask the difference in control law performance, and 
many flights must be made to achieve a valid comparison. In a digital simulation the noise in a pair 
of flights can be matched exactly with different control laws. Therefore, it is possible to apply 
statistical tests with sets of paired observations. The variation in performance for each pair due to 
turbulence variations is eliminated, which makes detection of performance differences easier. The 
variation in the turbulence for different pairs causes variations in the landing performance but only 
differences in performance due to different control laws are measured. The difference in the 
measured quantities for each pair of flights constitutes one observation, for example. 


y 



- ho2 


As usual, in testing the hypothesis that there is no difference in performance for the two control 
laws one applies Student’s t-test and calculates t (see, e.g., ref. 1 2) 


t 


= ^ ; 
/s^/n 




1 ^ 

n 




where n is the number of observations. For n=100 the critical regions are t> 1.9840 
and t<— 1.9840. When the resulting t is in one of the critical regions, we can say that with 
97.5 percent probability one of the systems performs better with respect to the variable y. Which 
system, of course^, is determined by the specific critical region, and the best estimate of the average 
improvement is y. 


AUTOMATIC CONTROL SYSTEM PERFORMANCE 


Comparison of Analog and Digital Control System Performance 

The basic digital control system was compared with an equivalent analog system. It is well 
known that the performance of a sample-data system becomes equal to that of an analog system as 
the sample interval goes to zero. Therefore, instead of developing a complete analog simulation 
including the aircraft equations of motion, the digital program was written in such a manner that 
the sample time interval could be made arbitrarily small. By reducing the interval from 0. 1 second 
(the interval used for most data runs) to 0.01 second, it was established that there were no 
significant differences in performance in calm or turbulent conditions. Few runs were needed, even 
in turbulence, since runs could be compared in matched pairs. The matched pairs were derived from 
solutions to the aircraft equations and stability augmentation equations, and from the turbulence 
samples at intervals of 0.01 second. Only the control computation cycle time was changed. 
Figure 15 shows a comparison of simulation flights. Since the turbulence is a function of altitude 
and time, slight variations in the flight path before the flare will change the turbulence encountered 
during the flare; and since the flare is sensitive to the exact form of the turbulence, the same 
turbulence is required for matched flights. The runs were, therefore, re-initialized at an altitude of 
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'>00 feet The flight paths were slightly different as would be expected from the difference in time 
intervals for control command changes. Statistically, however, there was no performance difference 
between the quasi-continuous and the digital system. 

Distance from glide-slope beam ground intercept, ft 



Figure 1 5.- Matched pair of flights in 20 ft/sec, headwind with turbulence. 


Performance of the Noncritical Control Modes 


In a comparison of the performance differences due to various modifications of the autopilot, 
the comparison must be made for each control mode separately. The altitude hold mode and the 
glide-slope capture mode were not altered after the first design since they performed satisfactorily 
under all conditions and they have the least influence on the final criterion for a successful landing. 
The flight paths in figure 16 for different glide-slope angles under no wind conditions illustrate the 
operation of these modes. Note the smooth glide-slope capture, with hardly any overshoots. The 
flight paths in the altitude hold mode, under no wind conditions, are simply straight horizontal 
lines. The glide-slope intercept altitudes were changed only to simplify the figure and are not 
intended to suggest deviation from the standard 1 500-foot approach. Figure 1 7 shows short sections 
of flights in turbulence for the altitude hold mode followed by the capture maneuver. Twenty flight 
paths are superimposed for each turbulent condition to show the widening of the flight envelope 
with increase in turbulence. Since low pitch rates are important for passenger comfort, the pitch 
command was limited to 2.5° per second in the altitude hold mode, but a rate limit was not 
imposed during the capture mode. Although not shown, simulation of the altitude hold mode with 
and without this limiter did not indicate any significant difference in flight paths. 
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Distance from ILS beam ground intercept, ft 


Figure 16.— Glide-slope capture for different glide-slope 
angles. 



Performance of the Critical Control Modes Under Steady Wind Conditions 

and Simple Shear 

Three configurations of the flare law were tested in steady winds with and without wind shear. 
In the first configuration (basic flare law), the predictive pitch command ramp (fig. 7) was applied 
as a function of time. The first modification was a pitch ramp based on altitude change that 
occurred in each sampling interval. The total pitch angle change was a predetermined constant. The 
second modification was the same as the first, except that vertical acceleration feedback was added. 
In the tests, the simulation was begun at an altitude of 1 500 feet about 9 miles from the runway in 
the automatic altitude hold mode. The wind velocity was varied between tests in increments of 
5 ft/sec, from a head wind of 20 ft/sec to a tail wind of 1 5 ft/sec. For the flights v^th a wind shear, 
a constant shear of 4 knots/100 ft was used. In all cases, with wind shear, the ground wind was zero. 

Steady winds had no significant effect on the performance of the altitude-hold, glide-slope 
capture, and glide-slope tracking modes. For the three flare modes, 1 iq, the sink rate at touchdown 
and Xq, the horizontal distance from the gliue slope at touchdown, were not significantly affected 
by steady wind, but Xq was reduced for the modified flare laws (figs. 1 8(a) and 1 8(b)). 

In the presence of wind shear, the landing performance of the three flare laws was similar to 
that without shear, but varied more with wind speed (figs. 19(a) and 19(b)). There was no 
significant effect on 1 iq, and the modified flare laws were superior for control of Xq, resulting in a 
350-foot shorter landing distance. 
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(a) Sink rate at touchdown vs. wind velocity; no 
wind shear. 



wind velocity, ft/sec 


(b) Horizontal distance from glide-slope beam 
center at touchdown vs. wind velocity; 
no wind shear. 


Figure 18.- Flare performance in the presence of steady wind. 



(a) Sink rate at touchdown vs. wind velocity; with 
wind shear of 4 knots/400 ft. 


(b) Horizontal distance from glide-slope beam 
center at touchdown vs. wind velocity 
with wind shear of 4 knots/ 100 ft. 


Figure 19.- Flare performance in the presence of wind shear. 
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Performance of the Critical Control Modes Under Turbulent Conditions 


The ultimate test of a system is its performance under the influence of natural disturbances. 
For the landing performance there is no one-to-one relation between the performance under steady 
conditions, as described in the last section, and performance under turbulent conditions. 

Glide-slope tracking— Table 1 shows the performance of the glide-slope tracking mode over a 
range of different wind conditions, where 

Ax average position error from the ideal position at 75 feet altitude on a 3° glide slope 
(x(actual) — x(ideal)) 

h^ average glide-slope tracking error (positive is above beam center) 

^h^ glide-slope tracking error standard deviation 

The simulation was initiated 45,000 feet from glide slope to ground intercept, and the above 
variables were calculated while the glide-slope tracking mode was engaged. 

TABLE 1 GLIDE-SLOPE TRACKING PERFORMANCE (100 RUNS PER CONDITION) 


Wind 


B 

% 

h=75ft 

h=7 5ft 

^Ax 

h=75ft 

he 


Percent 
outside 
FAA limit 

-20 

4 

10.97 


1.60 

52 


120 

-0.463 

13.1 

12.0 

-15 

3 

10.92 


1.20 

25 


110 

-.526 

10.03 

3.9 

-10 

2 

11.07 


.93 

7 


69 

-.962 

6.90 

.6 

-5 

1 

11.07 


,55 

8 


36 

-1.55 

3.72 

0 

5 

1 

11.64 


.47 

-61 


37 

.76 

3.56 

0 

10 

2 

1"" z3 


.91 

-130 


64 

.70 

6.6 

0 

15 

3 

x2.2Q 


1.20 

-197 


86 

1.29 

9.3 

4.9 


Examination of table 1 will disclose the following facts. The average glide-slope tracking error 
is slightly negative for head winds, which means that the airplane remains below the glide slope 
most of the time. The reverse is true for tail winds. The standard deviation of the tracking error 
increases with the standard deviation of the turbulence, as has already been shown qualitatively in 
figure 17. The average sink rate at 75 feet altitude is smallest for the weakest head wind and largest 
for the strongest tail wind because of the constant airspeed command of 230 ft/sec. The standard 
deviation of the sink rate at the flare window, is relatively small, and shows the expected trend 
of increase for turbulence. The percentage of flight time that the flight path was outside the limit 
specified in reference 1 1 is also shown. An on-board calculation of this error could be used to 
measure the severity of the turbulence, which would determine whether or not to continue with the 
automatic approach. The average position error Ax at 75 feet altitude follows the trend predicted 
from the average altitude error. It is, however, somewhat larger than would be predicted from the 
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altitude error along the glide slope alone. Here it must be remembered that the glide-slope tracking 
control begins to switch over to a sink rate control law below 200 feet altitude, and at 100 feet 
altitude no attempt is made to correct for beam tracking errors. Again, the standard deviation of the 
position errors at the flare window is largest for the greatest turbulence. 

The accuracy of achieving the flare window with the glide-slope extension was measured for a 
20 ft/sec head wind in turbulence. Without glide-slope extension the standard deviation of the rate 
of descent at 75 feet increased from 1.6 to 3.1 ft/sec, and the position error standard deviation 
increased from 1 20 to 260 feet. 

Flare mode- As discussed earlier, the flare maneuver is the most critical one for successful 
touchdowns in turbulence. For this reason various modifications to the flare control mode were 
studied. Table 2 summarizes the results for the more important changes investigated, each row 
summarizing a set of 100 touchdowns; the turbulence samples were matched for all flights, so that 
each set of flights could be compared to another set by means of the t-test described earlier. The 
table is subdivided into five groups according to the modifications investigated. This organization 
requires repetition of the data for some of the experiments, see column 2. 

A detailed explanation of the conditions for each experiment will be presented after some 
preliminary remarks about the significance of the data in some selected columns. Columns 1 to 5 
present experimental conditions and columns 6 to 16 present experimental results. The small tables 
at the bottom of the figure present statistical test results of comparing pairs of runs for sink rates, 
touchdown distances, and pitch at touchdown. The results are listed in table 3 with their standard 
deviations. 

The third parameter of importance in landings, the pitch angle at touchdown, Qq, and its 
standard deviation, will be discussed only briefly since, in all cases, the pitch angles were in the safe 
range to prevent a tail strike of the airplane. It is interesting that the differences in average pitch 
angle, although small, are statistically significant for many of the experiments. An examination of 
the raw data for each pair of flights showed that sign reversals in the term 0o, - ®02 were less 
frequent than sign reversals for hoj ^ ho 2 and Xq, — Xqj- For the vertical touchdown velocity our 
greatest interest is in the tail of the distribution, since it includes the cases of hard landings, which 
must be prevented. Therefore, a touchdown velocity frequency table, columns 14-16, is included 
for each experiment. 

The rows of table 2 present only a summary of the simulation results. The set of runs labeled 
Experiment 1 was the final “best system” selected after investigating various control law 
modifications. It will be noted that Experiment 1 is compared with other experiments for all except 
the first comparison, which was made to determine a control law configuration that was 
subsequently incorporated into the final system. Each subsequent comparison, tests 2-5, shows the 
effects of changing the control law from its configuration for Experiment 1 . In all cases the sink 
rate increases for a modification although sometimes not to a statistically significant degree. It 
would be generally desirable to have a short touchdown distance to leave adequate roll out distance 
on the runway. Low sink rate at touchdown and short touchdown distance are not naturally 
compatible; but since sink rate is considered to be more important, the best configuration chosen 
was a compromise weighted more heavily toward low sink rates. 
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TABLE 2.- COMPARISON OF CONTROL LAW MODIFICATIONS WHEN FLYING IN SEVERE TURBULENCE (a^, = 7.0 FT/SEC) 
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Statistically significant difference at the 0.995 significance level 
Statistically significaait difference at the 0.975 significance level 
No statistically significant difference 




The capability of selecting flare altitude and other flare parameters from the sink rate of the 
aircraft (Exp. 3) rather than having standard flare parameters (Exp. 8) was the first control 
improvement that was investigated. The average sink rate at touchdown was decreased by 
0.52 ft/sec and, as the frequency table shows, the variance was reduced by having fewer high sink 
rate touchdowns at the cost of adding distance to touchdown. It should be noted, however, that 
reducing the variation in touchdown distance resulted in better performance overall. For all further 
tests the sink rate dependent flare feature of Experiment 3 was incorporated in the final system. 

The second control modification studied was providing smaller integrator gain for pitch down 
than for pitch up; see column 4. The effect of this modification with otherwise the “best” flare law 
(Exp. I) is shown in section 2. The trend is similar to that for the first control law modification 
discussed. The average sink rate decreases and the average touchdown distance increases. 

The third control law modification studied was the effect of vertical acceleration feedback 
(column 3), A comparison of the sink rates for the best vertical acceleration feedback, nj\ = 0.9, and 
no acceleration feedback, nj^ = 1 .0, showed no statistically significant difference. Moreover, any 
difference would be of little practical interest, since the minimum is not very pronounced. An 
increase in vertical acceleration feedback (nj^ changing from 1.0 to 0.85) does decrease the average 
touchdown distance, while only slightly increasing its standard deviation. As a compromise, 
n^ = 0.9 was chosen. 

Pitch command rate limiting is required for passenger comfort. Since the flare maneuver 
requires greater pitch command activity, a pitch command rate limit of 5° per second was chosen, 
which is twice the limit of the pitch command rate for the glide-slope tracking mode. As can be seen 
from table 2, section 4, it was desirable to increase the rate limit, column 5. When the rate limit was 
completely removed (Exp. 6), the only statistically significant difference that occurred was an 
increase in touchdown distance. 

There is an inherent danger in the piecemeal improvement of a complicated system. If one 
optimizes one parameter and introduces additional control parameters, there is no assurance that 
the first parameter’s optimum adjustment will be optimum for the new system. Therefore, a final 
experiment was made (shown also in table 2, section 5), which proved that indeed there was a small 
advantage in computing the flare altitude, rather than preselecting it from the nominal glide-slope 
angle and approach speed, when otherwise the various improvements were kept. The improvement 
due to this portion of the flare law is not so clear cut as it was when it was first added to the analog 
equivalent system. 

A large sample of data for the 30 ft/sec head wind was taken in seven sections of 100 flights 
each and is represented in table 3. The purpose was to examine the variability between data runs 
with identical statistical turbulence, and to study the cause of high sink rate landings. The listing for 
each experiment of 100 runs includes a frequency table of the sink rates at touchdown. The last 
four lines of the table present the average of the seven runs, their standard deviation, the standard 
deviation divided by the average of each column, to show the relative variation of each column, and 
the probability distribution for the sink rate. In the touchdown sink rate frequency table, 
the ojyi data indicate that the center of the distribution is pretty well defined. However, the tails 
of the high sink rate distribution are not sufficiently well defined even with 700 individual flight 
records. For the design of actual flight systems, this lack of information would be serious, since we 
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are mainly interested in the number of hard landings that we can expect. It will also be noticed that 
the standard deviations of gust velocity aug = (2ai3g./100)'^^ averaged over each of the 100 flights 

are lower than would be expected from the gust model Uug = 0.2ug = 6 ft/sec (oug. is the variance 

of the gust about the mean for each particular run). This difference comes from the fact that for the 
given spectrum shape the average variance of a short sample is smaller than the variance of an 
infinite sample (see, e.g., ref. 1 2, fig. 3). 


TABLE 3.- DATA SUMMARY FOR SEVEN EXPERIMENTS OF 1 00 RUNS EACH IN SEVERE 
TURBULENCE (30 FT/SEC HEAD WIND, Uu = 6 FT/SEC) 


Experiment 

number 


ho 

no 

wind 

bo 

^^0 

^0 

no 

wind 

^0 

*0 


No. of touchdowns within indicatet 

sink rate limits, ft/s 

ec 


0-1 

1-2 

2-3 

3-4 

4-5 

5-6 

6-7 

7-8 

8-9 

9- 

10 

10- 

11 

1 

4.54 

2. 

02 

3.09 

1.3 

15 

60 

1666 
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4 

21 

23 

25 

20 

5 

1 

1 

0 
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2 

4.46 



2.84 

1.0 



1516 
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4 

18 

36 

29 

12 

0 

0 

1 

0 





3 

4,50 



3.12 

1.4 
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580 

5 

14 

33 

27 

12 

5 

2 

0 

2 





4 

4.61 



2.91 

1.2 



1660 
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4 

21 

31 

27 

12 

1 

4 

0 

0 





5 

4.43 



2.86 

1.2 



1558 
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3 

23 

34 

23 

13 

2 

2 

0 

0 





6 

4.68 



3.10 ' 

1.3 



1521 
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4 

20 

30 

22 

12 

10 

1 

1 
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7 

4.62 



2.89 

1,3 



1620 

1 
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5 

23 

28 

26 

10 

5 

2 

1 

0 





W 

4.55 



2.97 

1.24 



1592 

566.6 

4.14 

20 

30.7 

25.6 

13.0 

4.0 

1,71 

0.57 

0.258 
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0 

0 

0.0917 



0.124 

0.127 



61.8 ' 

54.2 

0.69 

3.1 

4,3 

2.4 

3.2 

3.3 

1.2 

0.53 

0.756 

{ 

3 

0 

Cj/u 

0.02 



0.042 

0.102 



0.038 

0.096 

0.16 

0.15 

0.14 

0.095 

0.247 
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0.0414 

0.2 
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0.04 

0.017: 

0.0057 
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We primarily want to know the major cause of hard landings. In the simulation it is easy to 
separate the effects of the vertical and horizontal gust components. One simply compares three sets 
of flights: one with both gust components present, one with vertical gust removed, and one with 
horizontal gust removed. As shown in table 4, the sink rate at touchdown was significantly smaller 
for the vertical gusts alone. Also, the dispersion of the sink rate at touchdown was decreased by a 
factor of 2. This indicates that for the assumed turbulence model horizontal gusts have a major 
effect on the sink rate. The horizontal gust patterns are shown in figure 20 for 1 7 flights with high 

TABLE 4.- EFFECT OF GUST COMPONENTS ON LANDING PERFORMANCE 
[25 flights each at 20 ft/sec head wind] 


Gust components 
present 

h 

% 

^o 

*^0 

Both components 

-2.4 

1.2 

1015 

580 

Horizontal only 

-2.5 

.99 

857 

460 

Vertical only 

-1.6 

.44 

901 

310 
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Figure 20.- Horizontal gust patterns that caused hard and soft landings; Ug^^ = 30 ft/sec, o = 6 ft/sec 





terminal sink rates and 1 7 flights with low terminal sink rates selected randomly from the 700 flight 
samples. It was pointed out earlier in the section pertaining to flights under steady wind conditions 
that the absolute magnitude of the wind is of little importance in the landing performance. For that 
reason, the horizontal reference line, in each graph, was drawn to intersect the ordinate at the wind 
velocity at touchdown. A vertical arrow indicates the time of the flare step command in pitch. In 
most hard landings first an increase and then a reduction of head wind occurs, meaning a sudden 
loss of lift. This also is indicated by the average wind velocity at touchdown of the 17 hard landings 
which IS 8.5 feet below the average wind velocity of 30 ft/sec. In contrast, the average wind at 
touchdown IS dose to the average wind velocity for the soft landings. Also, at touchdown, the wind 
velocity varies less and the large low frequency components are not present. 


Touchdown Performance of the Landing System 


The method of calculating the overall probability distribution of sink rate at touchdown is 
described in reference 13. Besides needing the probability distribution of ho at a given turbulence 
level, one also needs the probability distribution of the turbulence magnitude, provided there is 
turbulence. The only probability distribution of a^g available for landings was given in 
reference 13. A corrected version of this distribution, which had unit area as required for a 
probability distribution, was provided by Hawker, Siddeley Aviation, Ltd., through private 
communication (see fig. 21). This probability distribution was divided into eight intervals, and runs 
were made with Oug at the center of each interval. The probability distribution for ho for a given 
gust variance was weighted with the probability that this variance occurs, and the overall 
distribution of exceeding lio was calculated and is shown in figure 22 (curve 1) for the complete 
range for turbulent conditions shown in figure 21. For curve 2, landings in turbulence 
above Oug - 5 ft/sec are excluded. The results for both curves are similar to those given in 
reference 13, which includes data of actual landings. TTie dashed portions of the curves are 
extrapolations of the available data to a probability of 10 ®. ITiis value is a commonly used design 
goal for automatic landing systems. 



Figure 21 .- Horizontal gust distribution; all weather 
conditions, wind limit 20 knots. 
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Probability of excedance 



Figure 22 — Probability of exceeding a rate of descent 
(a straight line on this graph would indicate a 
normal distribution.) 


Clearly, if the results actually reflect 
reality, automatic landings with the system 
studied should not be attempted in more than 
moderate turbulence for landing gear damage 
probability to be less than 10 ®. This would 
mean that a certain percentage of all automatic 
landings should not be attempted because of 
turbulent conditions. 

Performance of the System With Reduced 
Sampling Rates 

It is of interest to know the minimum 
computation cycle rate, which does not 
compromise performance of the system. Table 5 
shows the performance of the final version of 
the system for different sample time intervals. 
For the interpretation of these data one must 
remember that the aircraft equations of motion 
and of stability augmentation were still solved at 
the rate of 1 0 times per second, and only the 
control equations were solved less often. Also, 
the parameters for the control equations were 
recomputed and adjusted for the specific sample 
time being tested (e.g., the coefficients of the 
difference equations are a function of the 
sampling time). 


Table 5(a) summarizes the glide-slope tracking performance changes with sampling time^^e 
j... ri Vi 'jt 7^ feet altitude (see footnote in table 5), arc closer to th 

Hare initiation conditions x and h than for the larger intervals, as seen 

A,», n,h,.pa.h e.„r standard 

deviation for tracking the glide slope, ah^, increases with increasing sampling time. 

Table 5(b) shows that for touchdown, as intervals become longer, there is a gradual decline m 
perfo™L ™fa^plane .ends to land harder and farther down 

pitch angle lends to be smaller. All three parameters have rncreased vanat A smk freq 
table is included to show the increase in hard touchdowns when 
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TABLE 5.- PERFORMANCE DEGRADATION DUE TO INCREASE IN SAMPLE TIME INTERVAL 

[vxg = -30 ft/sec] 

(a) Flight conditions before flare at h = 75 ft 


t 

7 1 

h 


_ 1 
X 




0.1 

-10.7 

1.4 

-1463 

130 

-0.98 

12.3 

.2 

-10.4 

1.9 

-1462 

140 

-1.0 

13.2 

.4 

-10.7 

2.2 

-1472 

150 

-1.4 

14.1 

.6 

-10.5 

2.2 

-1483 

140 

-0.8 

21.8 

.8 

-10.99 

2.3 

-1475 

160 

-2.1 

17.7 


' Ground track distance to ground intercept of glide slope 
and rate of descent at h = 75 ft, x = 1432.4 ft, and 
^nominal “ '10-4 ft/sec. 


(b) Touchdown 







No. 

of touchdowns within indicated 
sink rate limits (ft/sec) 

t 

ho 

no 

wind 

ho 


^0 

no 

wind 

^o 


e 

^0 

0-1 

1-2 

2-3 

3-4 

4-5 

5-6 

6-7 

7-8 

8-9 

0.1 

2. 

02 

2.77 

1.3 

15 

60 

937 

580 

0.0894 

0.0088 

5 ' 

27 

35 

15 

10 

6 

2 

0 

0 

.2 



2.78 

1.2 



1034 

620 

.0879 

.0084 

6 

26 

27 

26 

10 

3 

2 

0 

0 

.4 



12.90 

1.3 



1073 

630 

.0876 

.0088 

3 

20 

38 

21 

10 

3 

4 

1 

0 

.6 



3.13 

1.4 



1101 

660 

.0865 

.0091 

4 

21 

23 

27 

14 

7 

2 

2 

0 

.8 



3.73 

1.9 



1344 

710 

.0617 

.022 

5 

7 

28 

24 

15 

9 

4 

5 

2 
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Performance of the System With Randomly Missing Computation Cycles 


Randomly missing a few percent of the control computation cycles had little effect on the 
la,,d,no ptfo L" e A random number generator at each sample interval was used m dec.dmg 
wliclher or norto skip the control computations. An example with 5 percent mtssms samples is 
given in table 6. The effect is relatively small for the following two reasons. First, during le ™sse 
computation cycles, the control commands are simply frozen to the last commande va ue, an 
therefore not far from the correct value, and second, the system is a closed-loop feedbac sys em m 
which errors built up during a missed cycle are corrected by subsequent 

the extreme case for an average of 50 percent missing computation cycles the performance declines 
only slightly as evidenced by the increased touchdown distance dispersion. However the table for 
"nk rafe at’ Lchdown shows that hard landings will occur if the sampling rate is reduced 
randomly. This indicates that for the flare, samples should not be skipped. 


TABLE 6.- PERFORMANCE DEGRADATION DUE TO MISSING COMPUTATION CYCLES 


[t = 0.1 sec, 20 ft/sec head wind] 


(a) Flight conditions before flare at h - 75 ft 


Percent 

missing 

samples 

li 


X 


Se 

“e 

0 

10.7 

1.4 

1463 

130 

0.99 

12.3 

5 

11.06 

1.2 

1462 

110 

-0.899 

9.55 

50 

11.21 

1.8 

1462 

130 

-1.47 

12.9 


(b) Touchdown 



No. of touchdowns wii 
sink rate limits 

thin 

(ft/ 

indie 

sec) 

ated 


Percent 

missing 

samples 

li 

°h 

X 

Ox 

0-1 

1-2 

2-3 

3-4 

4-5 

5-6 

6-7 

7-8 

8-9 

0 

2.77 1 

1.3 

937 

580 

5 

27 

35 

15 

10 

6 

2 

0 

0 

5 

2.49 

1.0 

1375 

440 

8 

28 

36 

19 

6 

3 

0 

0 

0 

50 

2.76 

1.3 

1228 

880 

4 

27 

38 

18 

4 

6 

1 

1 

1 
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CONCLUSIONS 


The study has been concerned with automatic pitch axis control for the approach and landing 
of a large jet transport aircraft by means of a general purpose digital computer. The study 
considered the effects of wind shears and turbulence and assumed error free ILS data and 
knowledge of altitude, vertical velocity, vertical acceleration, pitch attitude, and pitch rate. 


The following conclusions were reached from this simulation study. 

1 . The chosen method of representing analog filters by difference equations was efficient and 
adequate. The frequency with which the control laws must be solved without loss of performance 
was relatively low, once every 0.2 second. (The required sampling rate for data acquisition and 
filtering may be appreciably higher.) 

2. Several flare laws implemented in this investigation decreased the variability of the landing 
performance in the face of gust disturbances. 

3. When the digital computer has other functions besides control (navigation, fuel 
management, measurement data smoothing and filtering, etc.), it is possible to skip randomly up to 
50 percent of the computation cycles without affecting the flight-path control, except the flare, to 
any great extent. 

4. Sink rate at touchdown is critically affected by transient wind shear starting during the 
flare. Therefore, an accurate gust gradient model of the turbulence near the ground is needed to 
predict the overall touchdown performance of automatic landing systems. Such a completely 
accurate gust model does not exist at this time. 

5. The Monte Carlo method of the study of hard landings presents a computer time problem. 
Since records of turbulence that cause hard landings are distinctively different from those that do 
not, it would be possible by means of a computer program to identify those simulated or measured 
turbulence records that may cause hard landings. If the probability of occurrence of such severe 
turbulence samples were determined as well, much computer time would be saved. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, June 19, 1970 
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APPENDIX A 


SIMULATION OF THE STABILIZED AIRPLANE 


In modem analog autopilots, control commands generated by the various autopilot control 
modes are not applied directly to the actuators. Instead, the control modes provide commands to a 
stabilized aircraft. In the case of the pitch axis control, pitch commands and velocity commands are 
supplied to the stabilization loop which in turn generates elevator and throttle position commands. 
In conventionally designed systems, the important state variables are fed back in a sequence of 
minor and major loops. Each successive feedback serves a different major function and involves 
different types of measurement and accuracy problems. This is in contrast to optimum control 
approaches where all state variables are fed back to meet a global performance criterion (see ref. 14, 
appendix E, for an optimal control approach to the Hare law design). The conventional approach 
has the advantage of simplicity and economy, since the airplane is already stabilized for all control 
modes. In this section, the aircraft, aircraft stabilization, and the resulting aircraft performance will 
be discussed. The system to be analyzed is shown in figure 1 . 

Both the airplane and the automatic landing system are simulated on the same digital 
computer. With the digital integration technique used, it was found that the airplane equations of 
motion had to be solved at least once each 0.1 second in aircraft time to result in a stable 
performance equivalent to an analog simulation. The Fortran IV simulation program was written so 
that the equations for analog system components were solved at the same rate as the aircraft 
equations of motion. The automatic control laws could be solved more slowly at any submultiple 
rate of the aircraft equations of motion. The solution rate of the aircraft equations of motion was 
adjusted so that any slow solution rate for the control equations could be simulated. The program 
was written either to run in real time or as fast as the computer could solve the equations. Fast runs 
at a computer cycle rate of 0.01 second per each 0.1 second aircraft time for both aircraft equations 
of motion and control equations allowed a 4-minute landing operation to be simulated in 
24 seconds. 


AIRPLANE EQUATIONS 


Equations of Motion 

The equations of motion of a jet transport in its landing configuration are represented in 
figure 1 by the block labelled aerodynamics, which shows the input and output variables. The 
stability derivatives and other characteristics are given in table 7. The airplane was considered to be 
representative of the current class of turbojet transport. Also, performance data were available from 
the analog simulation described in reference 2, which allowed comparison with the digital 
simulation for program verification. 

The equations of motion are simplified versions of those given in reference 2 modified to 
include the effect of wind. The modifications, which include removal of the lateral degrees of 
freedom and the effects of the landing gear contacting the ground, are straightforward and will not 
be discussed in detail. 
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Lift 


Cl(“) 

(ACl)ge 

0.19 + 5.3 a 
0.23 + 0.55 6f 
0.302 
7.68 

0.063(Cl^)£q£ 

Drag 1 

CoCa) 

0.03 + 0.08824 a + 1.8182 

(ACd)5^ 

(0.126 + 0.473 a)6f 

(^Co)gE 

(-0.02 - 0.332 a)e^^ 

be 

Pitch 

Cino 

0.09 


-1.062 

‘^”6e 

-0.923 

‘^”6f 

-0.103 

Cm • 

-4.01 

a 

Cmq 

-12.3 

(ACm)GE 

- 0 . 066 (CL„)cgg 

CACni)LG 

-0.01 

Ground effect factor 

^GE 

a972e"^^^^ 

h 

wheel weight 

Air density 

P 

2.377x10-3 slug/ft2 

Geometric parameters 

Area, S 

2758 ft2 

Chord, c 

22.16 ft 

Thrust offset, d 

4.0 ft 

Thrust angle, iy 

3.15'’ 

Weight parameters 

Weight, w 

200,000 lb 

Inertia, lyy 

3.9x106 slug-ft2 





The effects of wind are introduced into 
the equations as follows. First, local 
horizontal and vertical 
components Vxg and v^g of the wind are 

generated. These components are then 
transformed into two 
components Ug and wg in the trajectory axis 
system. This is a system whose x-axis (xj) is 

chosen to coincide with the tangent to the 
flight path (see fig. 23). (Ref. 15 states that 
flight simulator manufacturers often prefer 
this system of axes for exceptional stability.) 
As a second-order approximation, it is 
assumed that the wind acts on the aircraft as 
if it were a segment of the x-axis (ref. 16). 
Therefore, there is a contribution from the 
wind to the pitch rate as well as to the angle of attack and dynamic pressure. The effective pitch 
rate is given by equation (A6) where cog^ and cogy are the zj- axis components of gust at the 

wing and tail, respectively, and 80 feet is the distance between the aerodynamic centers of the wing 
and tail. In the simulation the tail gust was derived from the stored samples of the wing gust. The 
angle of attack due to the wind is Og of equation (A4), where Wg is the zj axis component of 
wind at the center of mass of the aircraft. The dynamic pressure, which affects the lift and drag, is 
calculated from the airspeed (eq. (A3)). Since stability derivatives are given for the wind axis system 
(x-axis coincided with the airspeed velocity vector), the lift and drag forces are transformed to 
components Fx and Fz in the trajectory axis system described above. 

The aircraft rate equations were numerically integrated using the fourth-order 
Adams-Bashford algorithm. With this algorithm a solution rate of t = 0.1 second was adequate for 
the basic aircraft equations of motion. For faster solution rates no significant differences in aircraft 
response were detected, while at slower rates the numerical solution became unstable. 

The equations of motion (Al) to (A16) are given in the order of the sequence, of calculations 
in the simulation program. The updating of the control inputs (6e, 5f, and T), which is a function of 
the automatic system and the servos, is not shown. The following assumptions are made in the 
aircraft equations of motion: 

1 . The airframe is assumed to be a rigid body with all motion restrained to the vertical plane. 

2. Earth is assumed to be flat, nonrotating, and fixed in inertial space. 

3. The mass of the airplane is assumed to be constant during each flight. 

4. Tire x-z plane is a plane of symmetry in the aircraft. 

5. Small angle approximations apply. 





Figure 23. - Vector diagram for equations of motion. 
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6. TTie air flow is assumed to be quasi-steady. 


7. 


Initial night conditions are for steady-state flight. Thus the problem can 
instantly m trimmed flight if the aircraft state variables are specified and the 
subroutines are initiated automatically. 


be started 
integration 


As shown in figure 23, four orthogonal axes systems are used in the same vertical plane: 

earth-fixed axis system; z points along the gravity vector; origin at the 
intersection of the glide slope and ground 




x,z 


xb^Zb 


trajectory axis system; reference frame with origin at aircraft center of mass 
and the x-direction alined with the tangent to the flight path 


wind axis system; aerodynamic frame of reference with the 
center of mass and the x-direction alined with the relative 


origin at the aircraft 
wind vector 


body axis system; x-axis fixed to the aircraft principal 
of the fuselage; origin at the aircraft center of mass 


axis along the centerline 


Figure 23 also identifies all other 
following set of definitions and equations. 


vector quantities and angles, which will appear in the 


Wind components in the earth axis system (average 
turbulence subroutine: 



wind + gust input) calculated in a 


CAl) 


Wind components along and perpendicular to the inertial velocity vector 

cos Y - Vz sin Y S Vx - Vz Y 
Wg = vxg sin Y + vzg cos Y = Vx^Y + vz^ 

Magnitude of the air velocity relative to the center of mass of the airplane: 

I V I = / (vg - Ug) + Wg2 
Angle between relative wind and inertial velocity vector: 



oig = sin' 


■1 


V 


(A4) 
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Aerodynamic angle of attack: 

a = n - «g 

Aerodynamic equivalent pitch rate due to gusts (ref. 16, p. 322). 

3Wg ^ ~ ^gT 

" 9x ' 80 

where 80 feet is the distance from wing to tail 


Aerodynamic pitch rate: 


^ " S 


Aerodynamic angle-of-attack rate: 

a = q - Y - q 

Pitch acceleration (inertial): 

pScV^ 


g 


q = 


2Iyy 


-ho 


Cm ot + Cm;,„'5e + 


•"66 


(ACm) 


(A5) 


(A6) 


(A7) 


(A8) 


pSc 
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YY 


V I t c„.i] 


dT 

Iyy 


CA9) 


Lift-to-mass ratio and drag-to-mass ratio in the wind axis system: 

k = ^ V^jcLCa) + CACL)^f + 


(AlO) 


D^^V2FcD(a) + (ACD)g£ + CACD)eEl ) 

m 2m L 

Aerodynamic force-to-mass ratios transformed through Og to the trajectory axis system. 


m 





\ 


> 


(All) 
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Acceleration components in the trajectory axis system: 


■ T Fx ~ T 

Vc = — cos Y - g sin Y - — = gY - 

m * mm* 


Fx 

m 


(A12) 


T 

azT = g cos Y - - sin (n + - -jj" = g 


‘T o --- ■ m 
Derivative of the flight-path angle: 


™ Cn + ix) - 


m 


Y = - 


Ve 


Integration of the inertial rates (fourth-order Adams-Bashford (see appendix B)): 

Y = J^y dt ; Vg = dt ; q = J" q dt ; Q = J' q dt 

Angle between the body axis and the tlight-path vector: 

n = 0 - Y 

Rate and position in the earth-fixed reference frame: 

fi = VgY 

h = y*h dt 


= J'w cos Y dt = J'' 


V dt 


m 


(A13) 

CA14) 

(AIS) 


(A16) 


This completes a computation cycle. The computation continues with the control equations 
and returns to equation (A1 ). 


Engine and Engine Servos 

Although the simulated reference jet transport represents a four-engine airplane, in the 
longitudinal simulation only, one equivalent engine needs to be simulated. The engine performance 
from throttle adjustment to change in thrust is represented as a first-order lag with a time constant 
of 1.25 sec. The engine servo is also assumed to be linear and it is represented as a second-order 
system with a natural frequency of 10 rad/sec, a damping ratio of 0.7, and a steady-state gain of 
1.0. This overall system from thrust command to response is represented as one difference equation 
in this simulation (see eq. (7) table 8). Advantages of this method are discussed in appendix B. 


Elevator, Elevator Servo, and Flaps 

The elevator and elevator servo combination is described by equation (7) table 8 also. In this 
case the elevator time constant is 0.2 second, natural frequency is 10 rad/sec, the damping ratio is 
0.7, and the steady-state gain is 0.8. 
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The flaps are represented as a strictly rate limited system. Upon a change in flap command, 
they advance at a rate of 0.035 rad/sec until the new command value is reached. 


AUTOMATIC FLIGHT-PATH STABILIZATION 

To control the flight path of an airplane automatically, it would be desirable to control the 
flight-path angle 7 directly. However, there is no output control variable that directly controls 7 . 
Also, for the basic airplane, for small deviations from horizontal flight, the steady-state flight-path 
angle is proportional to thrust for a constant pitch attitude, and the steady-state velocity decreases 
linearly with increase in pitch for constant thrust. Moreover, for short-period response, pitch does 
control the flight-path angle and pitch and pitch rate can be measured with gyros. To make short 
and long term responses agree, the aircraft is stabilized in the following manner. The airspeed of the 
aircraft is kept nearly constant by means of an autothrottle, and a pitch control loop maintains the 

flight path. 


Pitch Control Loop 

The pitch axis of the aircraft is controlled through an inner feedback loop to follow the pitch 
commands from the autopilot (fig. 1 ). The command signal is differenced with a pitch attitude 
signal for position feedback, and a filtered pitch rate signal for damping. The pitch rate filter, called 
a washout filter, responds quickly to a change in pitch rate, and does not respond to constant rates. 
The filter has been included for realism, since it would have been needed for banked turns which 
represent a constant pitch rate. The filter step response is k exp(-t/r), where k - 0.5 and r - 4 sec. 
Since there is no integration in the forward path of the pitch control loop, the error signal at the 
summing junction is not driven to zero for a given pitch command. Hence, for steady state, 
e = k0c where k is less than unity. This is alright in the automatic system, since the inner loop 
forms part of the forward path of the outer loop. Tlie outer loop does contain integration to null 
fliglit-path errors. In the simulation, the pitch control loop equations are solved at the same rate as 
the airplane equations of motion, since in an actual system this loop might remain an analog loop. 
The response of the airplane will be discussed following the description of the autothrottle. 


Autothrottle 

The autothrottle maintains constant airspeed by generating a thrust command signal to drive 
the throttle servo (fig. 24). The command signal is derived from three sources, the airspeed error, 
the longitudinal acceleration, and the pitch command. These three signals are further processed as 
follows. To provide an augmented airspeed error signal, the measured airspeed is summed with the 
airspeed reference and then summed with the output of a fore-aft accelerometer. The signal is then 
filtered in a blender filter (ref. 4). In effect, in turbulence the measured airspeed is heavily filtered 
so that the thrust command does not respond to gust inputs, whereas without turbulence the 
airspeed is not filtered. 
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Typical values: =4.0 sec 

ki ^ 0.097 sec/ft 
k| = 0.0048 sec/ft 
T 2 = 30.0 sec 
tq =2.0 sec 
k2 = 2.5 


Figure 24,— Autothrotlle. 


The response of the system described so far is not fast enough to maintain a constant airspeed 
during pitch maneuvers such as those encountered while tracking an ILS glide-slope beam. To 
compensate for this slow response, a predictive thrust command is derived from processing pitch 
command changes through a wash-out filter in sequence with a low pass filter. The augmented 
airspeed error signal, integral of airspeed error signal, and predictive thrust command are summed to 
form the thrust command signal. To prevent heavy throttle activity in turbulent conditions a rate 
limiter is required plus a position limit to prevent the thrust from being set dangerously low. 

In the simulation a steady-state thrust command will result in an equal engine thrust. In an 
actual autothrottic system only changes in throttle settings are commanded. Thrust magnitude is 
not commanded since actual thrust is not directly measured. However, the simulation is adequate, 
because changes in thrust settings are commanded only from the velocity errors and not from thrust 
errors. 
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APPENDIX B 


SIMULATION OF LINEAR SYSTEMS WITH DIFFERENCE EQUATIONS 


The traditional method of simulating linear systems by analog computer is to solve the 
differential equations of the system in terms of the highest derivative and integrate. However, in 
digital simulation, when the integration is carried out numerically, it can be a lengthy process that 
may be unstable for large sampling intervals. Simulating a linear system with the chosen difference 
equation is more efficient than using numerical integration techniques and is numerically stable 
independent of the sampling interval (ref. 17). 


The difference equation for one of the linear filters used in this simulation will be derived to 
provide a concrete example of the procedure. The transfer function of the wash-out filter used in 
the pitch stabilization loop is 


G(s) 


x(s) s + a 



The first step is to find the time response of the filter to a unit step input. 


yCs) 


ks ^ _ k 
s + a s s + a 


y(t) = ke-^^ 

Next, find the z transform of y(t). 



Factor the z transform of a unit step function out of the expression for y(z). 

z (unit step) = j- 


y(z) 


= k 


z - 1 


z - 


z 


z - 


1 


Since the response to a general input is wanted, replace the step input by a general input x(z). The 
actual input is approximated by a stair-step function that can be described by a linear combination 
of unit step functions and is equal to the general input at the sampling points. 


y(z) 


= k 


z - 1 


z 


e 


-aT 


x(z) 
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If wc cross multiply, 


(z - e"^^)y(z) = k(z - l)x(z) 

divide by z, 



transform to the discrete time domain, remembering that (l/z)yj^ = 

Yj\ ~ ^ yn-i “ ^ 

and rearranging, the final form of the difference equation is obtained: 

yn = e'^^Xn-1 + - Xn_i) 

At the sampling points, the difference equation is the exact solution for the response of the 
equivalent analog system to a stair-step function input. Therefore, the solution does not become 
numerically unstable for large sampling intervals. However, the output of the system will be realistic 
only if the sampling is done frequently enough to result in a stair-step waveform which is a good 
approximation of the original continuous input waveform. 

For comparison, the equations for the above filter using numerical integration (fourth-order 
Adams-Bashford) are: 

" ^n-1 ^ ~ " ^Xn-s) 

Xn = -ayn + kxn 

The computation efficiency of the difference equation is obviously better than that of the 
numerical integration technique. Also, this numerical integration method becomes numerically 
unstable for large sampling intervals, as confirmed during the investigation. 

Difference equations were used to simulate the aircraft control and automatic stabilization 
systems, to replace the analog autopilot with a digital version, and to generate the turbulent wind 
components. 

The aircraft control system simulator consists of a throttle servo with engine lag and an 
elevator servo with an elevator surface lag. Both of these servos have the same general form, a 
second-order transfer function with 0,7 critical damping in series with a first-order low pass filter. 
The form of the total transfer function and the corresponding difference equation are shown in 
section 7 of table 8. 

The washout filter in the stability augmentation system discussed above is simulated by the 
difference equation of section 1 in table 8. 
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All of the linear control laws in the autopilot were mechanized by difference equations to 
simulate a digital autopilot. The various control transfer functions are listed in sections 8 
tlirough 1 0 in table 8. 

The generation of the wind turbulence components requires the filtering of wide-band noise 
for spectrum shaping. Filters were simulated by the difference equations in sections 2 and 3 in 
table 8. 
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TABLE 8.- DIFFERENCE EQUATIONS USED IN THIS SIMULATION 




TABLE 8.- DIFFERENCE EQUATIONS USED IN THIS SIMULATION - Continued 
[Plus others of interest in control simulations] 



Era 


-AiAyAe + A7A5 + A3A7A4 



TABLE 8.- DIFFERENCE EQUATIONS USED IN THIS SIMULATION - Concluded 
[Plus others of interest in control simulations] 





APPENDIX C 


FILTER GAIN CONSTANT ADJUSTMENT FOR THE DIGITAL 
TURBULENCE SIMULATION 


For the turbulence simulation the variance of the gust must be generated correctly From th^ 
discussion of the difference equations in appendix B, the overall digital generation of the turbulence 
noise sequence can be compared to a continuous system with sample and hold circuits as shown in 
figure 25 The output of this circuit, expressed as a sequence of numbers, is statistically identica 
with the number sequence derived in the digital program /^erefore, the know edge can be used^ o 
calculate the ratio of output to input variance for both difference equations (eq ^( ) U ^ 

The input noise sequence to the filter can be thought of as a random wave (see fig. 26). Its power 

spectrum is (see ref. 9 and fig. 27). 2 

. . , _ ^x T j~ sin(a)T/2)1 
“ TT (,lT/2 


Gaussian 

random 

noise 



Outputs at 
the sampling 
times are 
identica! 


Figure 25.- Comparison of equivalent digital and analog turbulence generation. 
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Figure 26.— Random wave. 



Figure 27. Normalized power spectra, fI>(o) - 1 . 
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For example, it can be shown that at a sampling interval of t = 0.1 sec the input spectrum is flat to 
I dB up to 2.5 Hz. This is sufficiently high, since a transport aircraft approximated as a rigid body 
does not respond to such high frequencies, and since the gust spectrum «i>y(a;) is a low frequency 
spectrum (see fig. 27). Higher frequencies will be attenuated as a ’result of the nonwhiteness of the 
input spectrum. 


Tlie ratio of output variance to input variance is calculated from 


Oy 
^ a. 


f 


C^) I f (^) dco 


and 


/■ 


- / <t>x(w)d(i) 


From equation (1 6) of the text F(co) will be chosen so that 

1F(0))|^ = $y(0)) 

Using (Cl ) through (C4) gives the desired ratio 


(C2) 

(C3) 


(C4) 




2 


oo 

f (Ox^T/tt) j[sin((0T/2)]/(o)T/2) j%y(co)da 

00 

J* I ] / (wT/2) I d(jo 


(ox^T/^) J j[sin(wT/2)]/(o)T/2)p$y((o)dco 


CCS) 


with the approximation {[sin(cjT/2)]^o?T/2)}^ - 1 in the range where 4>y(a;) has significant energy 


2 


= CT/tt) J" 3>y(co)do) = (T/iT)ay^ 


CC6) 


An approach to exact equations for the output variance is to calculate it directly by expressing 
the difference equation involving input and output variables in terms of the input variable only. 

n 

Yn = / ^ CC7) 


1=0 


Since axj is the same for all xj, and the xj are independent, the variance of will be 




, . 2 


CCS) 


1=0 
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To find the aj one writes the initial terms of the difference equation for the first few inputs 
and then deduces the functional relationship between a, and aj+i . The correctness of this solution 
is then checked by induction. Now the sum of the squares of these terms can be split into geometric 
series and derivatives of geometric series, which are all summable. 

This method results in the following equation for 


2 

Cfu 


a 


iruu ^ 


e-auT)^ 

-2auT 
e “ 


(C9) 


Similarly, the development for begins by defining C ^e and rewriting (19), where D, 

and D 2 are defined as previously. 

yn = 2Cyn-i “ CVn-2 + D2^n-2 ^^10) 

With this, after some algebra, 






Dj£2. 


(1-C2)' 


+ Di2 


_1 

1 - C2 


2 


(Cii) 


For aT « 1, equations (C9) and (Cl 1) reduce to (C6). 
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